#### "Inclusionary regimes, party institutionalization and redistribution under authoritarianism" ####
# authors: "Lars Pelke"
# date: 2020-06-16
# written under "R version 3.6.1 (2019-07-05)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

#libraries
library(tidyverse)
library(ggpubr)
library(texreg)
library(dotwhisker)
library(broom.mixed)
library(countrycode)

library(lme4)
library(sjPlot)
library(sjmisc)
library(lfe)
library(sjPlot)
library(stargazer)
library(performance)
library(texreg)
library(ggeffects)
library(scales)


# set working directory
# please use the working directory, where you stored the zip-file. 

#### Load Data ####
vdem <- readRDS("data/vdem_merged.rds") 
vdem_spaw <- readRDS("data/vdem_merged_spaw.rds") 

vdem <- vdem %>%
  drop_na(v2x_regime)

vdem_spaw <- vdem_spaw %>%
  drop_na(v2x_regime)

#### Generate Sample ####

vdem <- vdem %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_spaw <- vdem_spaw %>%
  filter(v2x_regime == 0 | v2x_regime == 1)

vdem_rel <- vdem %>%
  drop_na(rel_red)

vdem_spaw <- vdem_spaw %>%
  drop_na(universalism_all)

vdem_rel <- distinct(vdem_rel, country_id, year, .keep_all= TRUE) # Distinct observations
vdem_spaw <- distinct(vdem_spaw, country_id, year, .keep_all= TRUE) # Distinct observations

vdem_rel <- vdem_rel %>% # delete those countries with less than 5 country-year observations to deal with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_spaw <- vdem_spaw %>% # delete those countries with less than 5 country-year observations to deal with singular fit 
  dplyr::group_by(country_id) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(country_id)

vdem_rel <- vdem_rel %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_spaw <- vdem_spaw %>%
  drop_na(pol_incl, v2xps_party, v2x_regime, year)

vdem_rel2 <- vdem_rel %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_spaw2 <- vdem_spaw %>%
  drop_na(e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, manufacturing_percent_ipol)

vdem_rel <- vdem_rel %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, rel_red, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_rel2 <- vdem_rel2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, rel_red, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

vdem_spaw <- vdem_spaw %>%
  select(country_id, country_name, year, pol_incl, v2xps_party, v2x_regime, universalism_all, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type)

vdem_spaw2 <- vdem_spaw2 %>%
  select(country_id, country_name, year, pol_incl, v2x_regime, 
         v2xps_party, universalism_all, e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, 
         v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, auto_regime_type, 
         urban_percent, manufacturing_percent_ipol)

#### Computing group-mean and de-meaned variables ####

vdem_rel <- sjmisc::de_mean(vdem_rel, pol_incl, v2xps_party, v2x_regime,
                            v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                            grp = country_id)

vdem_rel2 <- sjmisc::de_mean(vdem_rel2, pol_incl, v2xps_party, v2x_regime,
                             e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

vdem_spaw <- sjmisc::de_mean(vdem_spaw,pol_incl, v2xps_party, v2x_regime,
                            v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                            grp = country_id)

vdem_spaw2 <- sjmisc::de_mean(vdem_spaw2, pol_incl, v2xps_party, v2x_regime,
                             e_migdppcln, e_wb_pop_ln, v2exl_legitideolcr_1, urban_percent, manufacturing_percent_ipol, 
                             v2psorgs, v2psprbrch, v2psprlnks, v2psplats, v2pscohesv, 
                             grp = country_id)

#### Descriptive statistics, Appendix A1 ##

cols1 <- c("rel_red", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols2 <- c("rel_red", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols3 <- c("universalism_all", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm") 

cols4 <- c("universalism_all", "pol_incl_dm", "pol_incl_gm", 
           "v2xps_party_dm", "v2xps_party_gm", "v2x_regime_dm", "v2x_regime_gm", 
           "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm", "e_wb_pop_ln_gm",
           "v2exl_legitideolcr_1_dm", "v2exl_legitideolcr_1_gm", "urban_percent_dm", "urban_percent_gm", 
           "manufacturing_percent_ipol_dm", "manufacturing_percent_ipol_gm",
           "v2psorgs_dm", "v2psorgs_gm", "v2psprbrch_dm", "v2psprbrch_gm", 
           "v2psprlnks_dm", "v2psprlnks_gm", "v2psplats_dm", "v2psplats_gm","v2pscohesv_dm", "v2pscohesv_gm")

# Model 1 #
stargazer(
  title="Summary Statistics for Redistribution Dataset", 
  as.data.frame(vdem_rel[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Rel. Redistribution", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 2 #
stargazer(
  title="Summary Statistics for Universalism Dataset", 
  as.data.frame(vdem_spaw[, cols3]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Universalism", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", "Party Institutionalization (b)", 
                      "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 3 and 5 #
stargazer(
  title="Summary Statistics for Redistribution Dataset", 
  as.data.frame(vdem_rel2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Rel. Redistribution", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)

# Model 4 and 6 #
stargazer(
  title="Summary Statistics for Universalism Dataset", 
  as.data.frame(vdem_spaw2[, cols4]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"), 
  covariate.labels= c("Universalism", "Pol Inclusion (w)", " Pol Inclusion (b)",
                      "Party Institutionalization (w)", 
                      "Party Institutionalization (b)", "Electoral Autocracy (w)","Electoral Autocracy (b)", 
                      "GDP pc (w)", "GDP pc (b)", "Population log (w)", "Population log (b)", 
                      "Communist Ideology (w)","Communist Ideology (b)", 
                      "Urban Percentage (w)", "Urban Percentage (b)", "Manufacturing \ of GDP (w)", 
                      "Manufacturing % of GDP (b)",
                      "Party organizations (w)", "Party organizations", 
                      "Party branches (w)", "Party branches", "Party linkages (w)", "Party linkages", 
                      "Distinct party platforms (w)", "Distinct party platforms", 
                      "Legislative party cohesion (w)", "Legislative party cohesion")
)


## reference Level for Factor ##
vdem_rel$auto_regime_type  <- as.factor(vdem_rel$auto_regime_type)
vdem_rel2$auto_regime_type  <- as.factor(vdem_rel2$auto_regime_type)
vdem_rel$auto_regime_type = relevel(vdem_rel$auto_regime_type, ref=1)
vdem_rel2$auto_regime_type = relevel(vdem_rel2$auto_regime_type, ref=1)

vdem_spaw$auto_regime_type  <- as.factor(vdem_spaw$auto_regime_type)
vdem_spaw2$auto_regime_type  <- as.factor(vdem_spaw2$auto_regime_type)
vdem_spaw$auto_regime_type = relevel(vdem_spaw$auto_regime_type, ref=1)
vdem_spaw2$auto_regime_type = relevel(vdem_spaw2$auto_regime_type, ref=1)


#### Model Building step for step ####
## compare https://strengejacke.github.io/mixed-models-snippets/random-effects-within-between-effects-model.html 

## Growth Curve model

growth_model <- lmer(rel_red ~ year + pol_incl + v2xps_party +
                       auto_regime_type + (1 + year | country_id),
                     data = vdem_rel, 
                     REML = TRUE) # REML estimation for parameters
summary(growth_model)

## Complex REWB ##

m1<- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
            auto_regime_type + (1 + year | country_id) + (1 + pol_incl_dm + v2xps_party_dm | country_id),
          data = vdem_rel, 
          REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE

m1.1<- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
              auto_regime_type  + (1 + year | country_id) + (0 + pol_incl_dm + v2xps_party_dm | country_id),
            data = vdem_rel, # assume independence between random slopes and no coveriance
            REML = TRUE) # REML estimation for parameters
isSingular(m1.1, tol = 1e-05) # FALSE

## Simple REWB model ##

simpleREWB <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
                     auto_regime_type +  (1  | country_id),
           data = vdem_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(simpleREWB, tol = 1e-05) # FALSE

## Mundlak model ##

mundlak <- lmer(rel_red ~ year + pol_incl + pol_incl_gm + v2xps_party + v2xps_party_gm + 
                  auto_regime_type + (1 + year | country_id),
                data = vdem_rel, 
                REML = TRUE)
isSingular(mundlak, tol = 1e-05) # FALSE

#### Comparison of Models and check suitable model ####

tab_model(growth_model, m1, m1.1, simpleREWB, mundlak,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Classical RE", "Complex REWB", "Complex REWB (2)", "Simple REWB", "Mundlak"))

## REWB or simple RE-model ? ##
anova(growth_model, mundlak) # significant improvement of the Mundlak-model over the simple 
# RE-model, indicating that it makes sense to model within- and 
# between-subjects effects

## Simple or Complex REWB? ##

## REWB model or FE model ? ##

fixed <- felm(rel_red ~ year + pol_incl_dm + v2xps_party_dm +
                v2x_regime_dm | country_id,
              data = vdem_rel)

tab_model(fixed, growth_model, m1,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("FE-model with ID", "Classical RE", "Complex REWB"))

## simple REWB or Compley REWB? ##

anova(simpleREWB, m1) # no significant improvement of the complex REWB over the simple 
# REWB-model, indicating that it makes no sense to model random slopes and covariance

#### Residual diagnostics for complex REWB model ####
performance::check_model(simpleREWB)

#### Main Models ####

m1 <- lmer(rel_red ~  year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1, tol = 1e-05) # FALSE
performance::check_model(m1)


m2 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_spaw, 
           REML = TRUE) # REML estimation for parameters
isSingular(m2, tol = 1e-05) # FALSE
performance::check_model(m2)

m3 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m3, tol = 1e-05) # FALSE
performance::check_model(m3)


m4 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*auto_regime_type + v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_spaw2)
isSingular(m4, tol = 1e-05) # FALSE
performance::check_model(m4)


m5 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type +pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m5, tol = 1e-05) # FALSE
performance::check_model(m5)


m6 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2xps_party_dm + v2xps_party_gm + 
             auto_regime_type + pol_incl_dm*v2xps_party_dm*auto_regime_type +
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_spaw2)
isSingular(m6, tol = 1e-05) # FALSE
performance::check_model(m6)


tab_model(m1, m2, m3, m4, m5, m6, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1 Redistribution", "Model 2 Universalism", "Model 3", "Model 4", "Model 5", "Model 6"))

#### Regression Output ####

texreg(list(m1, m2, m3, m4, m5, m6), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:A11", 
       longtable = TRUE)

htmlreg(list(m1, m2, m3, m4, m5, m6), file = "Table1.doc", 
        inline.css = FALSE, doctype = TRUE, html.tag = TRUE, 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party Institutionalization",
                             "Party Institutionalization (b)",
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * PI",
                             "CM Autocracy  * PI", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "PI * Inclusiveness",
                             "HM Autocracy * PI * Inclusiveness",
                             "CM Autocracy * PI * Inclusiveness"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting relative redistribution",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "", # insert hmtl dagger here 
       label = "Table: A11")


##############################################################################################################
#### Margin Plots for Main Models ####
##############################################################################################################


#### Figure 1  based on Model 3 and 4 ####
#two way intercation effect plots 

margin_data <- ggpredict(m3, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  ylim(0,15) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m3, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  ylim(0,15) +
  labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m4, terms = c("pol_incl_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))

plot3 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  ylim(0,22) +
  labs(title = " ",
       x = "Pol Inclusion (within)",
       y = "Predicted universalism old-age pensions") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m4, terms = c("v2xps_party_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot4 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  ylim(0,22) +
    labs(title = " ",
       x = "Party Institutionalization (within)",
       y = "Predicted universalism old-age pensions") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


ggarrange(plot1, plot2, 
          labels = c("a", "b"), 
          ncol = 2, nrow = 1, 
          common.legend = TRUE)
ggsave("output/margins/Figure1_main.pdf", height = 12, width = 32, units= c("cm"))
ggsave("output/margins/Figure1_main.png", height = 12, width = 32, units= c("cm"), dpi = 1200)
ggsave("output/margins/Figure1_main_600dpi.png", height = 12, width = 32, units= c("cm"), dpi = 600)
ggsave("output/margins/Figure1_main_400dpi.png", height = 12, width = 32, units= c("cm"), dpi = 400)


ggarrange(plot3, plot4, 
          labels = c("a", "b"), 
          ncol = 2, nrow = 1, 
          common.legend = TRUE)
ggsave("output/margins/Figure1_SA.pdf", height = 12, width = 32, units= c("cm"))
ggsave("output/margins/Figure1_SA.png", height = 12, width = 32, units= c("cm"), dpi = 1200)


#### Figure 2 based on Model 5 and 6 #### 
# three way interaction effect plots 

sd(vdem_rel2$pol_incl_dm)

margin_data <- ggpredict(m5, terms = c("v2xps_party_dm", "pol_incl_dm[-0.37, 0.37]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.37", "0.37"),
                            labels = c("-1 SD", "+ 1 SD"))


plot1 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)



sd(vdem_spaw2$pol_incl_dm)

margin_data <- ggpredict(m6, terms = c("v2xps_party_dm", "pol_incl_dm[-0.43, 0.43]", "auto_regime_type"))

margin_data$facet <- factor(margin_data$facet, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Multiparty Autocracy"))
margin_data$group <- factor(margin_data$group, levels = c("-0.43", "0.43"),
                            labels = c("-1 SD", "+ 1 SD"))

plot2 <- plot(margin_data, 
              facet = TRUE) +
  labs(title = " ",
       color = "Political Inlusion (within)",
       x = "Party Institutionalization (within)",
       y = "Predicted universalism old-age pensions") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed") +
  scale_colour_grey(start = 0.1, end = 0.4) +
  scale_fill_grey(start = 0.1, end = 0.4)


ggarrange(plot1, 
          ncol = 1, nrow = 1, 
          common.legend = TRUE)
ggsave("Output/Margins/Figure2_main.pdf", height = 18, width = 32, units= c("cm"))
ggsave("Output/Margins/Figure2_main.png", height = 18, width = 32, units= c("cm"), dpi = 1200)
ggsave("Output/Margins/Figure2_main_600dpi.png", height = 18, width = 32, units= c("cm"), dpi = 600)
ggsave("Output/Margins/Figure2_main_400dpi.png", height = 18, width = 32, units= c("cm"), dpi = 400)


ggarrange(plot2, 
          ncol = 1, nrow = 1, 
          common.legend = TRUE)
ggsave("Output/Margins/Figure2_SA.pdf", height = 18, width = 32, units= c("cm"))
ggsave("Output/Margins/Figure2_SA.png", height = 18, width = 32, units= c("cm"), dpi = 1200)


#############################################################################################################
#### Main Models with Sub-Components of Party Institutionalization ####
#############################################################################################################


m1.1 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
             v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
             auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
             v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
             v2pscohesv_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_rel, 
           REML = TRUE) # REML estimation for parameters
isSingular(m1.1, tol = 1e-05) # FALSE
performance::check_model(m1.1)


m2.1 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
             v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
             auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
             v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
             v2pscohesv_dm*auto_regime_type + 
             (1  | country_id),
           data = vdem_spaw, 
           REML = TRUE) # REML estimation for parameters
isSingular(m2.1, tol = 1e-05) # FALSE
performance::check_model(m2.1)

m3.1 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m3.1, tol = 1e-05) # FALSE
performance::check_model(m3.1)


m4.1 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + pol_incl_dm*auto_regime_type + v2psorgs_dm*auto_regime_type + 
               v2psprbrch_dm*auto_regime_type + v2psprlnks_dm*auto_regime_type + v2psplats_dm*auto_regime_type + 
               v2pscohesv_dm*auto_regime_type + 
             e_migdppcln_dm + e_migdppcln_gm +
             e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
             urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_spaw2)
isSingular(m4.1, tol = 1e-05) # FALSE
performance::check_model(m4.1)


m5.1 <- lmer(rel_red ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
               v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
               v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_rel2)
isSingular(m5.1, tol = 1e-05) # FALSE
performance::check_model(m5.1)


m6.1 <- lmer(universalism_all ~ year + pol_incl_dm + pol_incl_gm + v2psorgs_dm + v2psorgs_gm + v2psprbrch_dm + v2psprbrch_gm +  
               v2psprlnks_dm + v2psprlnks_gm + v2psplats_dm + v2psplats_gm + v2pscohesv_dm + v2pscohesv_gm +
               auto_regime_type + v2psorgs_dm*auto_regime_type*pol_incl_dm + 
               v2psprbrch_dm*auto_regime_type*pol_incl_dm + v2psprlnks_dm*auto_regime_type*pol_incl_dm + v2psplats_dm*auto_regime_type*pol_incl_dm + 
               v2pscohesv_dm*auto_regime_type*pol_incl_dm + 
               e_migdppcln_dm + e_migdppcln_gm +
               e_wb_pop_ln_dm + e_wb_pop_ln_gm + v2exl_legitideolcr_1_dm + v2exl_legitideolcr_1_gm +
               urban_percent_dm + urban_percent_gm +  manufacturing_percent_ipol_dm +  manufacturing_percent_ipol_gm +
             (1  | country_id),
           data = vdem_spaw2)
isSingular(m6.1, tol = 1e-05) # FALSE
performance::check_model(m6.1)

tab_model(m1.1, m2.1, m3.1, m4.1, m5.1, m6.1, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1 Redistribution", "Model 2 Universalism", "Model 3", "Model 4", "Model 5", "Model 6"))

#### Regression Output ####

texreg(list(m1.1, m2.1, m3.1, m4.1, m5.1, m6.1), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Pol Inclusiveness",
                             "Pol Inclusiveness (b)",
                             "Party organizations", 
                             "Party organizations (b)", 
                             "Party Branches", 
                             "Party Branches (b)", 
                             "Party Linkages", 
                             "Party Linkages (b)", 
                             "Distinct Party Platforms", 
                             "Distinct Party Platforms (b)", 
                             "Legislative Party Cohesion", 
                             "Legislative Party Cohesion (b)", 
                             "Hegemonic Multiparty Autocracy",
                             "Competitive Multiparty Autocracy",
                             "HM Autocracy * Inclusiveness",
                             "CM Autocracy * Inclusiveness",
                             "HM Autocracy * Party organizations", 
                             "CM Autocracy * Party organizations", 
                             "HM Autocracy * Party Branches", 
                             "CM Autocracy * Party Branches", 
                             "HM Autocracy * Party Linkages", 
                             "CM Autocracy * Party Linkages", 
                             "HM Autocracy * Distinct Party Platforms", 
                             "CM Autocracy * Distinct Party Platforms", 
                             "HM Autocracy * Legislative Party Cohesion",
                             "CM Autocracy * Legislative Party Cohesion",
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Communist Ideology", 
                             "Communist Ideology (b)", 
                             "Urban Pop %",
                             "Urban Pop % (b)",
                             "Manufacturing Sector %",
                             "Manufacturing Sector % (b)",
                             "Inclusiveness * Party organizations",
                             "Inclusiveness * Party Branches",
                             "Inclusiveness * Party Linkages",
                             "Inclusiveness * Distinct Party Platforms",
                             "Inclusiveness * Legislative Party Cohesion",
                             "Inclusiveness * Party organizations * HM Autocracy",
                             "Inclusiveness * Party organizations * CM Autocracy",
                             "Inclusiveness * Party Branches * HM Autocracy",
                             "Inclusiveness * Party Branches * CM Autocracy",
                             "Inclusiveness * Party Linkages * HM Autocracy",
                             "Inclusiveness * Party Linkages * CM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * HM Autocracy",
                             "Inclusiveness * Distinct Party Platforms * CM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * HM Autocracy",
                             "Inclusiveness * Legislative Party Cohesion * CM Autocracy"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model predicting Redistribution / Universalism",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:A10", 
       longtable = TRUE)


#### Figure 3, Model 3 and 4 Interactions with different sub-components of Party Institutionalization ####

margin_data <- ggpredict(m3.1, terms = c("v2psorgs_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))


plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Organizations (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m3.1, terms = c("v2psprbrch_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Branches (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

# Plot 3
margin_data <- ggpredict(m3.1, terms = c("v2psprlnks_dm ", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot3 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Linkages (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

#Plot 4 
margin_data <- ggpredict(m3.1, terms = c("v2psplats_dm ", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot4 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Platforms (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

#Plot 5 

margin_data <- ggpredict(m3.1, terms = c("v2pscohesv_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot5 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Legislative Party Cohesion (within)",
       y = "Predicted relative redistribution") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


ggarrange(plot1, plot2, plot3, plot4, plot5,
          labels = c("a", "b", "c", "d", "e"), 
          ncol = 2, nrow = 3, 
          common.legend = TRUE)
ggsave("Output/Margins/Figure3_main.pdf", height = 25, width = 35, units= c("cm"))
ggsave("Output/Margins/Figure3_main.png", height = 25, width = 35, units= c("cm"), dpi = 1200)
ggsave("Output/Margins/Figure3_main_600dpi.png", height = 25, width = 35, units= c("cm"), dpi = 600)
ggsave("Output/Margins/Figure3_main_400dpi.png", height = 25, width = 35, units= c("cm"), dpi = 400)



#### Model 4.1 ####

margin_data <- ggpredict(m4.1, terms = c("v2psorgs_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))


plot1 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Organizations (within)",
       y = "Predicted Universalism") +
  geom_hline(yintercept=0, linetype="solid") +
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


margin_data <- ggpredict(m4.1, terms = c("v2psprbrch_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot2 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Branches (within)",
       y = "Predicted Universalism") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

# Plot 3
margin_data <- ggpredict(m4.1, terms = c("v2psprlnks_dm ", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot3 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Linkages (within)",
       y = "Predicted Universalism") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

#Plot 4 
margin_data <- ggpredict(m4.1, terms = c("v2psplats_dm ", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot4 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Party Platforms (within)",
       y = "Predicted Universalism") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")

#Plot 5 

margin_data <- ggpredict(m4.1, terms = c("v2pscohesv_dm", "auto_regime_type"))

margin_data$group <- factor(margin_data$group, levels = c("0", "1", "2"),
                            labels = c("Closed Autocracy", "Hegemonic Multiparty Autocracy", "Competitive Mulitparty Autocracy"))

plot5 <- plot(margin_data, 
              colors = "bw",
              facet = TRUE) +
  labs(title = " ",
       x = "Legislative Party Cohesion (within)",
       y = "Predicted Universalism") +
  geom_hline(yintercept=0, linetype="solid") + 
  geom_hline(yintercept=median(vdem_rel2$rel_red), linetype="dashed")


ggarrange(plot1, plot2, plot3, plot4, plot5,
          labels = c("a", "b", "c", "d", "e"), 
          ncol = 2, nrow = 3, 
          common.legend = TRUE)
ggsave("Output/Margins/Figure3_SA.pdf", height = 25, width = 35, units= c("cm"))
ggsave("Output/Margins/Figure3_SA.png", height = 25, width = 35, units= c("cm"), dpi = 1200)


#### Empirical Examples for introduction and empirical section ####


vdem_cases <- vdem_rel2


vdem_cases <- vdem_cases %>%
  select(country_name, year, pol_incl, pol_incl_dm, pol_incl_gm, v2xps_party, 
         v2xps_party_dm, v2xps_party_gm, auto_regime_type, rel_red) %>%
  filter(country_name %in% c("Vietnam", "Singapore", "Belarus")) 

vdem_cases %>%
  filter(country_name == "Vietnam") %>%
  summarize(Mean = sd(v2xps_party))


#### Model Fit with Empirical Cases: Vietnam, Singapore and Belarus ####

#### M3 ####

vdem_rel2 <- vdem_rel2 %>%
  drop_na(manufacturing_percent_ipol_dm)

vdem_rel2$fitted_values_m3 <- fitted(m3)

cases_df <- vdem_rel2 %>%
  dplyr::select(country_name, country_id, year, v2xps_party_dm, v2xps_party_gm, v2xps_party, pol_incl, pol_incl_dm, pol_incl_gm, 
         auto_regime_type, rel_red, fitted_values_m3) %>%
  filter(country_name %in% c("Vietnam", "Singapore", "Belarus")) %>%
  mutate(diff_pred = rel_red - fitted_values_m3)

cases_df %>%
  group_by(country_name) %>%
  summarize(mean = mean(diff_pred), 
            sd = sd(diff_pred))

cases_df_belarus <- cases_df %>%
  filter(country_name == "Belarus")

cases_df_singapore <- cases_df %>%
  filter(country_name == "Singapore")

cases_df_vietnam <- cases_df %>%
  filter(country_name == "Vietnam")

## Belarus ##

plot_belarus <- ggplot(data = cases_df_belarus, aes(x = year)) +
  geom_line(aes(y = rel_red, color = "Relative Redistribution"), size = 1.2) +
  geom_line(aes(y = fitted_values_m3, color = "Fitted Values"), size = 1.4) +
  geom_line(aes(y = v2xps_party*20, color = "Party Institutionalization"), linetype = "dotted", size = 1.2) +
  geom_line(aes(y = pol_incl*20, color = "Political Inclusion"), linetype = "twodash", size = 1.2) +
  scale_y_continuous("Redistribution", sec.axis = sec_axis(~./20, name = "Party Institutionalization /\n Polical Inclusion")) +
  theme_pubr() +
  labs(color = "", 
       title = "Belarus")

##Singapore##
plot_singapore <- ggplot(data = cases_df_singapore, aes(x = year)) +
  geom_line(aes(y = rel_red, color = "Relative Redistribution"), size = 1.2) +
  geom_line(aes(y = fitted_values_m3, color = "Fitted Values"), size = 1.4) +
  geom_line(aes(y = v2xps_party*20, color = "Party Institutionalization"), linetype = "dotted", size = 1.2) +
  geom_line(aes(y = pol_incl*20, color = "Political Inclusion"), linetype = "twodash", size = 1.2) +
  scale_y_continuous("Redistribution", sec.axis = sec_axis(~./20, name = "Party Institutionalization /\n Polical Inclusion")) +
  theme_pubr() +
  labs(color = "", 
       title = "Singapore")

## Vietnam ##
plot_vietnam <- ggplot(data = cases_df_vietnam, aes(x = year)) +
  geom_line(aes(y = rel_red, color = "Relative Redistribution"), size = 1.2) +
  geom_line(aes(y = fitted_values_m3, color = "Fitted Values"), size = 1.4) +
  geom_line(aes(y = v2xps_party*20, color = "Party Institutionalization"), linetype = "dotted", size = 1.2) +
  geom_line(aes(y = pol_incl*20, color = "Political Inclusion"), linetype = "twodash", size = 1.2) +
  scale_y_continuous("Redistribution", sec.axis = sec_axis(~./20, name = "Party Institutionalization /\n Polical Inclusion")) +
  theme_pubr() +
  labs(color = "", 
       title = "Vietnam")

ggarrange(plot_singapore, plot_belarus, plot_vietnam,
          ncol = 1, nrow = 3, 
          common.legend = TRUE)
ggsave("Output/A11.pdf", height = 30, width = 25, units= c("cm"))
ggsave("Output/A11.png", height = 30, width = 25, units= c("cm"), dpi = 1200)

